import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pyswarm import pso
from scipy import stats
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.model_selection import train_test_split
from tqdm import tqdm

np.random.seed(42)

plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号

# 载入实验数据
file_path = 'D:\\workapp-for-computer\\111code\\2024-9-21-huawei-bei\\code-c\\data\\附件一（训练集）.xlsx'
data = pd.read_excel(file_path, sheet_name='材料1')

# 筛选出只含有正弦波的行
data = data[data['励磁波形'] == '正弦波']

# 数据提取
temperatures = data.iloc[:, 0].values  # 温度
frequencies = data.iloc[:, 1].values  # 频率
losses = data.iloc[:, 2].values  # 磁芯损耗
B_p_values = data.iloc[:, 4:1028].values  # 磁通密度数据

# 选择一个磁通密度峰值
B_p = np.max(B_p_values, axis=1)

# 数据探索
print("原始数据描述统计:")
print(data.describe())

# 绘制箱线图
plt.figure(figsize=(10, 6))
sns.boxplot(data=data[['温度，oC', '频率，Hz', '磁芯损耗，w/m3']])
plt.title('原始数据箱线图')
plt.show()

# 使用 Z-Score 方法识别异常值
z_scores = stats.zscore(data[['温度，oC', '频率，Hz', '磁芯损耗，w/m3']])
abs_z_scores = np.abs(z_scores)
filtered_entries = (abs_z_scores < 3).all(axis=1)

# 过滤掉异常值
cleaned_data = data[filtered_entries]

# 再次绘制箱线图查看效果
plt.figure(figsize=(10, 6))
sns.boxplot(data=cleaned_data[['温度，oC', '频率，Hz', '磁芯损耗，w/m3']])
plt.title('Z_score处理后数据箱线图')
plt.show()

# 保存处理后的数据
output_file_path = 'D:\\workapp-for-computer\\111code\\2024-9-21-huawei-bei\\code-c\\data\\process_附件一.xlsx'
cleaned_data.to_excel(output_file_path, index=False, sheet_name='材料1')
print(f"处理后的数据已保存至 {output_file_path}")

# 重新载入处理后的数据
processed_data = pd.read_excel(output_file_path, sheet_name='材料1')

# 筛选出只含有正弦波的行
processed_data = processed_data[processed_data['励磁波形'] == '正弦波']

# 数据提取
temperatures = processed_data.iloc[:, 0].values  # 温度
frequencies = processed_data.iloc[:, 1].values  # 频率
losses = processed_data.iloc[:, 2].values  # 磁芯损耗
B_p_values = processed_data.iloc[:, 4:1028].values  # 磁通密度数据

# 选择一个磁通密度峰值
B_p = np.max(B_p_values, axis=1)

# 定义斯坦麦茨方程
def steinmetz_eq(params, B_p, f):
    k, a, b = params
    return k * (B_p ** a) * (f ** b)

# 定义新的温度修正方程##############################################################
def new_temp_adjusted_eq(params, B_p, f, T):
    k1, alpha_1, beta_1, gamma, eta = params
    tau_0 = 15  # 假设参考温度为25°C
    return k1 * f ** alpha_1 * B_p ** beta_1 * (1 + gamma * (T - tau_0) + eta * (T - tau_0) ** 2)

# 计算均方误差
def mse(y_true, y_pred):
    return mean_squared_error(y_true, y_pred)

# 计算平均绝对误差
def mae(y_true, y_pred):
    return mean_absolute_error(y_true, y_pred)

# 粒子群算法目标函数，计算斯坦麦茨方程的预测误差
def objective(params):
    losses_pred = steinmetz_eq(params, B_p_train, frequencies_train)
    return mean_absolute_error(losses_train, losses_pred)

# 粒子群算法目标函数，计算新的温度修正方程的预测误差
def objective_new_temp(params):
    losses_pred = new_temp_adjusted_eq(params, B_p_train, frequencies_train, temperatures_train)
    return mean_absolute_error(losses_train, losses_pred)

# 参数的上下界
lb = [0, 0, 0]  # k, a, b的下限
ub = [100, 5, 5]  # k, a, b的上限
lb_new_temp = [0, 0, 0, -0.1, -0.1]  # k1, alpha_1, beta_1, gamma, eta 的下限
ub_new_temp = [100, 5, 5, 0.1, 0.1]  # k1, alpha_1, beta_1, gamma, eta 的上限

# 按温度分组并划分训练集和测试集
temperature_groups = cleaned_data.groupby('温度，oC')

temperature_results = {}
temperature_mae_results = {}

total_mae_raw = 0
total_mse_raw = 0
total_mae_new = 0
total_mse_new = 0

for temperature, group in tqdm(temperature_groups, desc="按温度分组拟合斯坦麦茨方程"):
    # 划分训练集和测试集
    train_group, test_group = train_test_split(group, test_size=0.3, random_state=42)

    # 数据提取
    B_p_group = np.max(train_group.iloc[:, 4:1028].values, axis=1)
    f_group = train_group.iloc[:, 1].values
    l_group = train_group.iloc[:, 2].values
    T_group = train_group.iloc[:, 0].values

    # 使用粒子群算法拟合斯坦麦茨方程
    optimal_params, optimal_value = pso(
        lambda params: mean_absolute_error(l_group, steinmetz_eq(params, B_p_group, f_group)), lb, ub, maxiter=200,
        swarmsize=50)

    # 使用粒子群算法拟合新的温度修正方程
    optimal_params_new_temp, optimal_value_new_temp = pso(
        lambda params: mean_absolute_error(l_group, new_temp_adjusted_eq(params, B_p_group, f_group, T_group)),
        lb_new_temp, ub_new_temp, maxiter=5000, swarmsize=1000)

    k_pso, a_pso, b_pso = optimal_params
    mse_st = mse(l_group, steinmetz_eq(optimal_params, B_p_group, f_group))
    mae_st = mae(l_group, steinmetz_eq(optimal_params, B_p_group, f_group))

    k1_pso, alpha_1_pso, beta_1_pso, gamma_pso, eta_pso = optimal_params_new_temp
    mse_new_temp = mse(l_group, new_temp_adjusted_eq(optimal_params_new_temp, B_p_group, f_group, T_group))
    mae_new_temp = mae(l_group, new_temp_adjusted_eq(optimal_params_new_temp, B_p_group, f_group, T_group))

    temperature_results[temperature] = {
        'steinmetz_params': optimal_params,
        'new_temp_params': optimal_params_new_temp,
        'steinmetz_mae': mae_st,
        'steinmetz_mse': mse_st,
        'new_temp_mae': mae_new_temp,
        'new_temp_mse': mse_new_temp
    }

    temperature_mae_results[temperature] = {
        'steinmetz_mae': mae_st,
        'steinmetz_mse': mse_st,
        'new_temp_mae': mae_new_temp,
        'new_temp_mse': mse_new_temp
    }

    total_mae_raw += mae_st
    total_mse_raw += mse_st
    total_mae_new += mae_new_temp
    total_mse_new += mse_new_temp

# 打印不同温度下的拟合结果
for temperature, result in temperature_results.items():
    print(f"温度 {temperature}°C:")
    print(
        f"斯坦麦茨方程粒子群算法拟合的参数: k={result['steinmetz_params'][0]:.4f}, a={result['steinmetz_params'][1]:.4f}, b={result['steinmetz_params'][2]:.4f}")
    print(f"斯坦麦茨方程的平均绝对误差: {result['steinmetz_mae']:.4f}")
    print(f"斯坦麦茨方程的均方误差: {result['steinmetz_mse']:.4f}")
    print(
        f"新的温度修正方程粒子群算法拟合的参数: k1={result['new_temp_params'][0]:.4f}, alpha_1={result['new_temp_params'][1]:.4f}, beta_1={result['new_temp_params'][2]:.4f}, gamma={result['new_temp_params'][3]:.4f}, eta={result['new_temp_params'][4]:.4f}")
    print(f"新的温度修正方程的平均绝对误差: {result['new_temp_mae']:.4f}")
    print(f"新的温度修正方程的均方误差: {result['new_temp_mse']:.4f}\n")

print(f"原始方程总的mae: {total_mae_raw:.4f}")
print(f"原始方程总的mse: {total_mse_raw:.4f}")
print(f"修正方程总的mae: {total_mae_new:.4f}")
print(f"修正方程总的mse: {total_mse_new:.4f}\n")

# 使用粒子群算法拟合斯坦麦茨方程（整体数据）
train_data, test_data = train_test_split(processed_data, test_size=0.3, random_state=42)

# 训练集数据提取
temperatures_train = train_data.iloc[:, 0].values
frequencies_train = train_data.iloc[:, 1].values
losses_train = train_data.iloc[:, 2].values
B_p_train_values = train_data.iloc[:, 4:1028].values
B_p_train = np.max(B_p_train_values, axis=1)

# 测试集数据提取
temperatures_test = test_data.iloc[:, 0].values
frequencies_test = test_data.iloc[:, 1].values
losses_test = test_data.iloc[:, 2].values
B_p_test_values = test_data.iloc[:, 4:1028].values
B_p_test = np.max(B_p_test_values, axis=1)

# 拟合整体数据的斯坦麦茨方程
with tqdm(total=200, desc="斯坦麦茨方程拟合进度") as pbar:
    def objective_with_progress(params):
        loss = objective(params)
        pbar.update(1)
        return loss

    optimal_params, optimal_value = pso(objective_with_progress, lb, ub, maxiter=200, swarmsize=50)
k_pso, a_pso, b_pso = optimal_params
print(f"斯坦麦茨方程粒子群算法拟合的参数: k={k_pso:.4f}, a={a_pso:.4f}, b={b_pso:.4f}")
print(f"斯坦麦茨方程的平均绝对误差: {optimal_value:.4f}")
mse_st = mse(losses_train, steinmetz_eq(optimal_params, B_p_train, frequencies_train))
print(f"斯坦麦茨方程的均方误差: {mse_st:.4f}")

# 使用粒子群算法拟合新的温度修正方程
with tqdm(total=5000, desc="新的温度修正方程拟合进度") as pbar:
    def objective_new_temp_with_progress(params):
        loss = objective_new_temp(params)
        pbar.update(1)
        return loss

    optimal_params_new_temp, optimal_value_new_temp = pso(objective_new_temp_with_progress, lb_new_temp, ub_new_temp,
                                                          maxiter=5000, swarmsize=1000)
k1_pso, alpha_1_pso, beta_1_pso, gamma_pso, eta_pso = optimal_params_new_temp
print(
    f"新的温度修正方程粒子群算法拟合的参数: k1={k1_pso:.4f}, alpha_1={alpha_1_pso:.4f}, beta_1={beta_1_pso:.4f}, gamma={gamma_pso:.4f}, eta={eta_pso:.4f}")
print(f"新的温度修正方程的平均绝对误差: {optimal_value_new_temp:.4f}")
mse_new_temp = mse(losses_train, new_temp_adjusted_eq(optimal_params_new_temp, B_p_train, frequencies_train, temperatures_train))
print(f"新的温度修正方程的均方误差: {mse_new_temp:.4f}")

# 斯坦麦茨方程预测损耗
losses_st_train = steinmetz_eq((k_pso, a_pso, b_pso), B_p_train, frequencies_train)
losses_st_test = steinmetz_eq((k_pso, a_pso, b_pso), B_p_test, frequencies_test)

# 新的温度修正方程预测损耗
losses_new_temp_adj_train = new_temp_adjusted_eq((k1_pso, alpha_1_pso, beta_1_pso, gamma_pso, eta_pso), B_p_train, frequencies_train, temperatures_train)
losses_new_temp_adj_test = new_temp_adjusted_eq((k1_pso, alpha_1_pso, beta_1_pso, gamma_pso, eta_pso), B_p_test, frequencies_test, temperatures_test)

# 可视化结果 (散点图)
plt.figure(figsize=(12, 6))
# plt.scatter(range(len(losses_train)), losses_train, label='训练集实际损耗', marker='o', color='blue')
# plt.scatter(range(len(losses_st_train)), losses_st_train, label='训练集斯坦麦茨方程预测', marker='x', color='red')
# plt.scatter(range(len(losses_new_temp_adj_train)), losses_new_temp_adj_train, label='训练集新的温度修正方程预测', marker='^', color='green')

plt.scatter(range(len(losses_test)), losses_test, label='测试集实际损耗', marker='o', color='blue')
plt.scatter(range(len(losses_st_test)), losses_st_test, label='测试集斯坦麦茨方程预测', marker='x', color='red')
plt.scatter(range(len(losses_new_temp_adj_test)), losses_new_temp_adj_test, label='测试集新的温度修正方程预测', marker='^', color='green')

plt.xlabel('样本索引')
plt.ylabel('磁芯损耗 (瓦特每立方米)')
plt.title('磁芯损耗预测比较 (散点图)')
plt.legend()
plt.grid(False)
plt.show()



# 提取温度和对应的误差值
temperatures = list(temperature_mae_results.keys())
mae_st_values = [result['steinmetz_mae'] for result in temperature_mae_results.values()]
mse_st_values = [result['steinmetz_mse'] for result in temperature_mae_results.values()]
mae_new_values = [result['new_temp_mae'] for result in temperature_mae_results.values()]
mse_new_values = [result['new_temp_mse'] for result in temperature_mae_results.values()]

# 创建散点图
plt.figure(figsize=(12, 6))

# 原始方程 MSE 散点图
plt.scatter(temperatures, mse_st_values, color='red', label='原始方程均方误差 (MSE)', marker='o')

# 原始方程 MAE 散点图
plt.scatter(temperatures, mae_st_values, color='red', label='原始方程平均绝对误差 (MAE)', marker='^')

# 修正方程 MSE 散点图
plt.scatter(temperatures, mse_new_values, color='green', label='修正方程均方误差 (MSE)', marker='o')

# 修正方程 MAE 散点图
plt.scatter(temperatures, mae_new_values, color='green', label='修正方程平均绝对误差 (MAE)', marker='^')

plt.title('原始方程与修正方程的误差比较')
plt.xlabel('温度 (°C)')
plt.ylabel('误差值')
plt.xticks(rotation=45)  # 旋转x轴标签
plt.legend()
plt.grid(False)
plt.tight_layout()  # 调整子图间距
plt.show()
